Tutorial Part II — R goes spatial

Introduction

Let me start with a big thanks to Moritz Wenzler, Cedric Scherer, Pierre Gras and Juergen Niedballa who constantly helped to improve the course! I am deeply indepted to you!

The course is intended to be a quick introduction to using R as a GIS. It builds on the previous course, Day1_RIntro, and is written in such a way that you can understand handling spatial data when knowing the basics of handling data.frames, matrices and lists. There is no need to know tidyverse/dplyr coding style.

This course provides the basics for using the most important spatial data types — vector (points, lines, polygons, aka ESRI shapefiles) and raster data. In R, these are Spatial(Points/Lines/Polygons)DataFrame or RasterLayer Objects. Recently, sf-objects have been developed to handle vector data more easily. The latter can be treated like simple data.frames. Also the raster package is under development — called terra package - but the course does not cover this. The terra package is especially useful for remote sensing data.

In addition to different data types, the course provides sections on coordinate projections and the most important geospatial operations. Last but not least - we plot a lot to ease data visualisation and because it is fun.

The data we use stem from a longstanding project we run in Borneo. Please refer to our website www.ecological-dynamics-izw.com. The project involves species conservation and large scale landscape planning. To learn more about the data used in the course and the project, please refer to the following publications that are freely available:

Targeted conservation to safeguard a biodiversity hotspot from climate and land-cover change. MJ Struebig, et al. 2015. Current Biology 25 (3), 372-378
Anticipated climate and land‐cover changes reveal refuge areas for Borneo’s orang-utans. MJ Struebig, et al. 2015. Global Change Biology 21 (8), 2891-2904
The importance of correcting for sampling bias in MaxEnt species distribution models. S Kramer‐Schadt et al. 2013. Diversity and Distributions 19 (11), 1366-1379.
The Borneo carnivore database and the application of predictive distribution modelling. S Kramer-Schadt, et al.2016. Raffles Bulletin of Zoology, Supplement No. 33: 18–41.

Useful (web)sites

R news and tutorials

http://www.r-bloggers.com/

Quick introduction to spatial R

Spatial visualisation

check out the packages {tmap}, {ggmap}, {leaflet} and {cartography}
https://cran.r-project.org/web/packages/tmap/vignettes/tmap-getstarted.html
A whole book: https://clauswilke.com/dataviz/ with spatial data here: https://clauswilke.com/dataviz/geospatial-data.html

overview of spatial packages

R-cheatsheets — great for learning ‘vocabulary’

in your folder
or here: https://www.rstudio.com/resources/cheatsheets/

Infos about simple features in R ({sf} package)

https://r-spatial.github.io/sf/index.html
and nice intro: https://oliviergimenez.github.io/introspatialR/#1
geocomputation: https://geocompr.robinlovelace.net/spatial-operations.html

Use EPSG codes as unique and specific identifies of your coordinate reference systems instead of writing projection details. See http://spatialreference.org/ for EPSG codes or https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf
BioMove specials:

For analysis of plot-based biodiversity data: package {vegan}
For analysis of telemetry data: packages {adehabitatHR}, {adehabitatLT}, {move}, {recurse}, {momentuHMM}, moveHMM},{ctmm}`

Basics

How to install…

pkgs = c("here", "sp", "sf", "dismo", "raster", "GISTools",  "rgdal", 
         "maptools", "rgeos","rgl","rasterVis", "ellipse", "units", "tmap", 
         "viridis","dplyr", "gdistance","vegan", "patchwork",
         "psych", "iNEXT", "BAT", "ggplot2", "stars") # package names
# install.packages(pkgs)
# update.packages()

…and load packages

library(here) #for easy directory management

# soon not going to be supported any ore
library(sp) 
library(raster)
library(rgdal)
library(rgeos)
library(rgl)
library(rasterVis)

# the soon up-to-date spatial packages: 'terra' and 'stars'
# however, they still create problems when used together with 'raster' etc.
library(sf)
library(stars)
library(tmap)
library(dismo)
library(GISTools)
library(maptools)

library(viridis) # nice colour palettes
library(units)
library(ggplot2)
library(patchwork) # to combine plots

How to call help

# Information about a function or a package. If you do not understand some
# of the code below, always check the arguments and help of the function.
?dismo
# search for an item:
??SpatialPolygons
# show the instructions of a package:
vignette(package = "dismo") # vignette(package="sp")
# load the instruction of a package (embedded pdf):
vignette("over") #vignette("sdm")

sp::plot
showMethods(f = 'plot') # S4 type method
methods(generic.function = 'plot') # S3 type method

Assign the workspace

Set your root directory, e.g. the directory under which you store everything you need for this course, i.e. the data and the R-Scripts. You could name it ‘TUB_BD_Course’.

Create the subfolders relative to root-folder. Please adjust to your setup. A possible folder structure could be:

├── TUB_BD_Course                   # root folder 
│   └── _Data_Borneo                # the Borneo data
│       └── BaseMaps                # geo data, found in _Data_Borneo
│       └── BaseRecords             # animal observation data in _Data_Borneo
│   └── R                           # store here all your scripts, i.e.
│       └── script_course1.R
│       └── script_course2.R

Traditionally, the working directroy was hard coded, i.e. the full path was specified:

# adjust the working directory [wd]
## getwd()  # alternative fct
## 
## # Set your OWN root directory, the following is just an example of how to do it:
## root_wd <- setwd("C:/Users/kramer/PopDynIZW Dropbox/Steph Kramer/_GitHub")

This approach is error-prone and it complicates the cooperation with others as they have to update the working directory in every script they receive from you. (And if they don’t change it back to your directory, you have to update it again…)

Nowadays, the best approach is to work with R projects that are associated with own working directory, workspace, history, and source documents.

You can create a project in the RStudio IDE by using the Create Project command (available on the Projects menu and on the global toolbar). In our case, we want to create a project in an existing directory where there is already R code and data placed inside that said folder. Now, a hidden folder called .Rproj.user and most importantly a file called your_porject.Rproj shoudl appear in your folder. Every time you want to work on your project, you open theRstudio session by double-clicking that .Rproj file which ensures that the root working directory is correctly set.

When using R projects, the {here} is a handy helper:

Illustration here package by Allison Horst

The goal of the {here} package is to enable easy file referencing in project-oriented workflows. In contrast to using setwd(), which is fragile and dependent on the way you organize your files, here uses the top-level directory of a project to easily build paths to files.

The here() function from the {here} package searches for the .Rproj file and will allow you to construct relative paths easily within your project environment:

here::here()

## [1] "C:/Users/kramer/PopDynIZW Dropbox/Steph Kramer/_GitHub"

here::dr_here()

root_wd <- here()

Now, we assign the path to these folders relative to the root directory:

├── root_wd                    # TUB_BD_Course 
│   └── dbor_wd                # _Data_Borneo
│       └── maps_wd            # BaseMaps
│       └── recs_wd            # BaseRecords

#today, we start with Borneo
dbor_wd     <- here("_DataBorneo")
  maps_wd   <- paste0(dbor_wd, "/BaseMaps") # mind difference paste0() and paste()
  recs_wd   <- paste(dbor_wd, "/BaseRecords", sep = '')

Create an output folder where you can temporarily store your results (anything you plot, create, want to save,….) with dir.create.

#output_wd <- paste0(dirname(root_wd), "/", "output")
output_wd <- here("output")
if (!dir.exists(output_wd)) {dir.create(output_wd)}

Finally, your folder structure will look like this:

├── TUB_BD_Course                   # root folder 
│   └── _Data_Borneo                # the Borneo data
│       └── BaseMaps                # geo data, found in _Data_Borneo
│       └── BaseRecords             # animal observation data in _Data_Borneo
│   └── R                           # store here all your scripts, i.e.
│       └── script_course1.R
│       └── script_course2.R
│   └── output                      # for any results/ plots,...

Access folder content

You will find a lot of bioclimatic raster files (.asc) as well as shapefiles in the BaseMaps folder:

filenames <- list.files(path = maps_wd)
## run also:
# filenames
head(x = filenames)

## [1] "bio_asc_01.asc" "bio_asc_02.asc" "bio_asc_03.asc" "bio_asc_04.asc"
## [5] "bio_asc_05.asc" "bio_asc_06.asc"

Recap

useful functions:
head() : only shows the first 6 elements, i.e. if you have a large data frame, you can quickly check header and the first entries (analogously: tail() )
paste(x, sep='') : appends characters or numbers — great for working directories, loops across many maps or data.frames etc
data formats and access:
data.frame() : [row, column] or “dollar sign”, i.e. mydata[1,2]: first row, second column, or mydata$MyColumnName[1]
matrix() : see above
list() : [[element]], i.e. mylist[[1]]. Attention! A list element can consist of a data.frame, that you can then access like that: mylist[[1]][1,2]…
S4-objects (like spatial data) : access via @, i.e. RasterLayer@data@values

Data Types

Raster data

Data import

First, load a raster map. Please check the description in the folder called DataDescription_Borneo.doc. The file >bio_asc_01.asc is showing mean annual temperature * 10. Why? Because it can be stored as integer, not floating type, which saves a lot of storage and memory.

## # read the ascii file as raster format
ras_bio_asc_01 <- raster(x = paste0(maps_wd, "/bio_asc_01.asc"))

## or use here:
ras_bio_asc_01 <- raster(x = here("_DataBorneo", "BaseMaps", "bio_asc_01.asc"))

Bor_mat <- ras_bio_asc_01 # easy copying of whole maps

Now have a look at the object

ras_bio_asc_01

## class      : RasterLayer 
## dimensions : 1425, 1423, 2027775  (nrow, ncol, ncell)
## resolution : 0.008333334, 0.008333334  (x, y)
## extent     : 108.3341, 120.1925, -4.374013, 7.500988  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : bio_asc_01.asc 
## names      : bio_asc_01

What does the slot with extent tells you? Looks like decimal degrees, i.e. could be latitude and longitude. However, the slot crs (crs stands for coordinate reference system) is empty. That means, without knowing the reference system, the map cannot be plotted at the exact location on the globe. Thus, when you get a map, always make sure you also get the information of the coordinate reference system! If it is not yet set (because you import a raster map as a simple matrix or points from a file), this is the first thing you should do! Btw — you can access the information of an S4-class such as RasterLayer objects with @:

# str(ras_bio_asc_01) #check this
ras_bio_asc_01@crs    # crs = coordinate reference system. If empty, assign it!

## Coordinate Reference System:
## Deprecated Proj.4 representation: NA

Assign coord. ref. system

Important! Please read and repeat basic GIS-knowledge on coordinate reference systems (CRS), map projections, and geographic and projected coordinates.

If a map already has a defined CRS, then you cannot simply overwrite this! To transform the coordinates, e.g. angular coordinates of latitude/ longitude of WGS84, into a projected CRS with planar coordinates (e.g. equal area of Lambert Azimuthal), you need to transform it (see below chapter Raster Projection).

In the following, we assign the CRS only when no information is provided. This can be the case when you upload a map based on a simple data frame or matrix, i.e. when it is not a spatial object.
@crs gives the information of the object (here: raster layer) Bor_mat or ras_bio_asc_01, and CRS() is calling to a function that sets/defines the CRS of the object. Mind the spelling in small letters/ capitals!

crs = shows coordinate reference system of a spatial object
CRS() creates projection and sets arguments for CRS!

An easy way to specify the CRS is using the EPSG code, e.g. instead of CRS("+proj=longlat +datum=WGS84") use CRS("+init=epsg:4326").

See http://spatialreference.org/. The EPSG code for WGS84 is 4326.

ras_bio_asc_01@crs <- CRS(projargs = "+proj=longlat +datum=WGS84")# equal to
ras_bio_asc_01@crs <- CRS("+init=epsg:4326")
# mind: crs != CRS()

Projections can also be transferred from one object to another, but only when the CRS hasn’t yet been specified, i.e. when the @crs-slot is empty. In the following, load the topographic map (digital elevation model) of Borneo (Bor_dem.asc or bio_asc_24.asc) and assign it the same CRS as Bor_mat, because we know they are overlaying (you can check that with the raster extents):

# ras_bio_asc_24 = DEM = digital elevation model
ras_bio_asc_24 <- raster(x = paste0(maps_wd, "/bio_asc_24.asc")) 
# or use here:
# ras_bio_asc_24 <- raster(x = here("_DataBorneo", "BaseMaps", "bio_asc_24.asc")) 

ras_bio_asc_24@crs <- ras_bio_asc_01@crs # same as: CRS("+proj=longlat +datum=WGS84")

Plot for a quick look

There are different ways to plot a raster map. Below is a chapter on visualising raster data, here I just give a very quick overview to get a first impression of the data we work with.

The easy base plot:

# base plot
plot(x = ras_bio_asc_01)

Using the {tmap} package:

# static plot
tmap_mode(mode = "plot")
tm_shape(shp = ras_bio_asc_01) + tm_raster()

Very good for checking the correct location of the data. Click on the + | - | layer icons on the upper left and zoom in or out and change the background layers:

# interactive plot
tmap_mode(mode = "view")
tm_shape(shp = ras_bio_asc_01) + tm_raster()

Clip areas

We do this first to work with a small raster to save computation time. For this, we create a clip_extent based on the spatial coordinates. The command ‘crop’ then clips the raster map.

extent(x = ras_bio_asc_01) # same as extent@ras_bio_asc_01

## class      : Extent 
## xmin       : 108.3341 
## xmax       : 120.1925 
## ymin       : -4.374013 
## ymax       : 7.500988

clip_extent <- extent(117.2, 117.45, 5.4, 5.5)
ras_bio_asc_01_cr <- crop(x = ras_bio_asc_01, y = clip_extent)
plot(ras_bio_asc_01_cr, col = inferno(10))

Accessing RasterLayers

First, have a look at the internal data structure of the raster object:

ras_bio_asc_01_cr # str(ras_bio_asc_01_cr) # attributes(ras_bio_asc_01_cr) # class(ras_bio_asc_01_cr)

## class      : RasterLayer 
## dimensions : 12, 30, 360  (nrow, ncol, ncell)
## resolution : 0.008333334, 0.008333334  (x, y)
## extent     : 117.2008, 117.4508, 5.400988, 5.500988  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : bio_asc_01 
## values     : 261, 268  (min, max)

There are many ways to retrieve internal data and access single bits of information:

bbox(obj = ras_bio_asc_01_cr)

##           min        max
## s1 117.200790 117.450790
## s2   5.400988   5.500988

ras_bio_asc_01_cr@extent@xmin

## [1] 117.2008

# xmin(ras_bio_asc_01_cr) # xmin(extent(ras_bio_asc_01))
# ras_bio_asc_01_cr@extent # extent(ras_bio_asc_01_cr) # attr(ras_bio_asc_01_cr,'extent')
ras_bio_asc_01_cr@ncols        # ncol(ras_bio_asc_01_cr) !n.b. cols vs. col!

## [1] 30

head(ras_bio_asc_01_cr@data@values)  # values(ras_bio_asc_01_cr)

## [1] 265 265 265 266 267 267

ras_bio_asc_01_cr@crs

## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World."],
##         BBOX[-90,-180,90,180]]]

# crs(ras_bio_asc_01_cr) # projection(ras_bio_asc_01_cr) !=CRS(ras_bio_asc_01_cr)

# crs = coordinate reference system defined
# CRS creates projection and takes args for crs!
# e.g. CRS("+proj=longlat +datum=WGS84")

#Retrieve internal data cont.:
nrow(x = ras_bio_asc_01_cr) # ncol(ras_bio_asc_01_cr) # ncell(ras_bio_asc_01_cr)

## [1] 12

dim(x = ras_bio_asc_01_cr)  # 12 rows, 30 columns, 1 z-dimension

## [1] 12 30  1

res(x = ras_bio_asc_01_cr)  # resolution = cell size

## [1] 0.008333334 0.008333334

An important function is coordinates(), which returns the centroids of each raster cell:

head(x = coordinates(obj = ras_bio_asc_01_cr), n = 10) # command from package 'sp'

##              x        y
##  [1,] 117.2050 5.496821
##  [2,] 117.2133 5.496821
##  [3,] 117.2216 5.496821
##  [4,] 117.2300 5.496821
##  [5,] 117.2383 5.496821
##  [6,] 117.2466 5.496821
##  [7,] 117.2550 5.496821
##  [8,] 117.2633 5.496821
##  [9,] 117.2716 5.496821
## [10,] 117.2800 5.496821

Terrain computation

We will now work with the digital elevation model (DEM) ras_bio_asc_24 and create a hillshade (3D look) based on topography using slope and aspect: - simulates a 3D surface - computes shaded relief values for a raster surface

slope <- terrain(x = ras_bio_asc_24, opt = "slope",
                 unit = "radians", neighbors = 8)
aspect <- terrain(x = ras_bio_asc_24, opt = "aspect",
                  unit = "radians", neighbors = 8)
Bor_hs <- hillShade(slope, aspect, angle = 45,
                    direction = 270)

Let’s plot the hillshade:

plot(Bor_hs, col = grey(0:100/100), legend = FALSE)

Other very useful terrain calculations are cost surfaces, cost distances and least cost path, e.g. for corridor calculations. We will do that in another course.

Visualising rasters

Difference between plot and image

plot

plot keeps the proportion of the map/ cells (e.g. squares here).

plot(x = ras_bio_asc_01, col = viridis(256))

#plot(ras_bio_asc_01_cr)

However, if you plot several layers on top of each other, plot will not overlay them exactly, because the argument ‘add’ fails (example below). In that case, better use image, or any other plot function.

image (raster)

If you use image, the plotted proportions will be changed and skewed to the plot surface.

image(x = ras_bio_asc_01)

# image(ras_bio_asc_01_cr)

`{tamp}` (interactive plots)

You can also save the interactive map — check your output folder for the html file. Tip of the day: always save the epsg code of your crs at the end of a spatial file.

tmap_mode("view")
(i.m <- tm_shape(ras_bio_asc_01) +
  tm_raster(palette = "viridis",  title = "Mean annual temp"))

tmap_save(i.m, paste0(output_wd,"/BorneoMAT_4326.html"))
# or use here: 
# tmap_save(i.m, here("output","BorneoMAT_4326.html"))

You can also use this package to create high quality plots for your publications.

tmap_mode("plot")
(m <- tm_shape(ras_bio_asc_01) +
  tm_raster(palette = "viridis",  title = "Mean annual temp"))

tmap_save(m, paste0(output_wd,"/BorneoMap_4326.png"), 
          units = "mm", width = 70, height = 70, dpi = 150 )

Plot raster data with `{ggpolot2}`

For advanced R-users: Beautiful plots can also be made with {ggpolot2}; however, RasterLayer objects cannot directly be plotted, but have to be converted to a conventional data frame first, i.e. 2 steps are needed prior to plotting. The object SpatialPointsDataFrame is explained in more detail when working with vector data.

# First, convert RasterLayer object to a SpatialPointsDataFrame
ras_bio_asc_01_pts <- rasterToPoints(ras_bio_asc_01, spatial = TRUE)
# Then to a dataframe
df_bio_asc_01  <- data.frame(ras_bio_asc_01_pts)
class(df_bio_asc_01)

## [1] "data.frame"

head(df_bio_asc_01)

##   bio_asc_01        x        y optional
## 1        271 116.7383 7.030154     TRUE
## 2        271 116.7466 7.030154     TRUE
## 3        270 116.7300 7.021821     TRUE
## 4        270 116.7383 7.021821     TRUE
## 5        270 116.7466 7.021821     TRUE
## 6        270 116.7550 7.021821     TRUE

ggplot(data = df_bio_asc_01) +
  geom_raster(aes(x = x, y = y, fill = bio_asc_01)) + 
  coord_sf(crs = "+init=epsg:4326") + ## set correct projection
  scale_fill_viridis(name = "°C * 10") + ## set custom colors
  labs(x = "Longitude", y = "Latitude",
       title = "Mean annual temperature") +
  theme_minimal() ## set custom plot style

Alternatively, one can use the {stars} package. (The {stars} package is another package aimed at working with spatial data, namely spatiotemporal arrays such as raster and vector data cubes). We first turn opur raster into a stars object and can then use the geom_stars() function in combination with ggplot:

stars_bio_asc_01 <- st_as_stars(ras_bio_asc_01)
class(stars_bio_asc_01)

## [1] "stars"

ggplot() +
  geom_stars(data = stars_bio_asc_01) +
  scale_fill_viridis(name = "°C * 10", na.value = "transparent") + ## set custom colors
  labs(x = "Longitude", y = "Latitude",
       title = "Mean annual temperature") +
  theme_minimal() ## set custom plot style

Composite plots

Simple composite plot

Plot the hillshade (3D relief) and add temperature colours on top: alpha value gives semi-transparency.

image(Bor_hs, col = grey(0:100/100))
image(ras_bio_asc_24, col = terrain.colors(25, alpha = 0.3), add = TRUE)
points(coordinates(ras_bio_asc_01_cr), cex = 0.1, pch = '+')
plot(extent(ras_bio_asc_01_cr), add = TRUE, col = 'red')

Now open the image-plot with the zoom icon in a separate window and change the size. You will see that the different layers are still plotted on top of each other.

Now do the same with the base plot function:

plot(Bor_hs, col = grey(0:100/100), legend = FALSE)
plot(ras_bio_asc_24, col = terrain.colors(25, alpha = 0.3), add = TRUE)
points(coordinates(ras_bio_asc_01_cr), cex = 0.1, pch = '+')
plot(extent(ras_bio_asc_01_cr), add = TRUE, col = 'red')

The extent is plotted as a red box, and the centroid coordinates of the raster cells are plotted as points. Use the ‘zoom’ function to increase the plot. To do so, call ‘zoom’ and click twice on the map to span your rectangle. Check out the area around the red box:

# zoom()
# locator(2) # gives you the coordinates of — here 2 — locations interactively

Composite plot with `{tmap}`

In the following, we will create a SpatialPointDataFrame from the small clipped/ cropped raster with the {sf} package and plot it on top of the hillshade. Use the possibility of selecting/ unselecting different layers (layer icon on the left).

hillsh <- Bor_hs

# make sf object from coordinates -> you will learn that also later
ras_bio_asc_01_cr_sf <- st_as_sf(
  data.frame(coordinates(ras_bio_asc_01_cr)), # create dataframe of coordinates
  coords = c("x", "y"), # define columns for the coordinates
  crs = 4326, # define crs, 4326 is the EPSG code
  sf_column_name = "geometry" # sf needs a geometry column and you have to name it
)


# interactive plotras_bio_asc_24
tmap_mode(mode = "view")
tm_shape(shp = hillsh, raster.downsample = TRUE)  +  
    tm_raster(palette = "Greys") +
  tm_shape(shp = ras_bio_asc_24, raster.downsample = TRUE) + 
    tm_raster(palette = grDevices::topo.colors(20),alpha = 0.3) +
  tm_shape(shp = ras_bio_asc_01, raster.downsample = TRUE) + 
    tm_raster(palette = grDevices::rainbow(10), alpha = 0.3) +
  tm_shape(shp = ras_bio_asc_01_cr_sf) + 
  tm_dots(shape = 20, size = 0.01)

What do the shape-arguments mean, e.g. shape = 20? Have a look here at the symbols: https://ggplot2.tidyverse.org/articles/ggplot2-specs.html?q=shape#sec:shape-spec

Unfortunately, in {tamp} you can only use shape = 1 and shape = 20, only circles are plotted.

# static plot  
tmap_mode(mode = "plot")
tm_shape(shp = hillsh, raster.downsample = TRUE) + 
    tm_raster(palette = "Greys") +
  tm_shape(shp = ras_bio_asc_24, raster.downsample = TRUE) +
    tm_raster(palette = grDevices::terrain.colors(10), alpha = 0.3) +
  tm_shape(shp = ras_bio_asc_01_cr_sf) + 
  tm_dots(shape = 1, size = 0.05)

perspective: 3D plot

# Cool 3D plots with rgl library, e.g. 'rgl.surface'
persp(x = ras_bio_asc_24, xlab = "Easting", ylab = "Northing",
      zlab = "elevation", r = 5, d = 1.5, expand = 0.1,
      ticktype = "detailed")

Why is it so spiny??? Check units! This is still a geographic crs, not a projected one. Wait for the lecture for details.

Raster projection

Project the grid (the small one!!) to have all units in meters:

proj_moll = CRS("+proj=moll +lat_0=65 +lon_0=10") # Mollweide projection
Bor_dem_moll <-projectRaster(ras_bio_asc_01_cr, crs = proj_moll)
persp(Bor_dem_moll, xlab = "Easting", ylab = "Northing",
      zlab = "elevation", main = "Elevation model of Borneo",
      r = 1, d = 5.5, expand = 0.1, ticktype = "detailed")

Raster Stacks

A raster stack is a collection of RasterLayer objects with the same spatial extent and resolution, similar to a geodatabase. Go into the maps folder and check what is inside:

# gives names and full path of file
files.full <- list.files(path = maps_wd, pattern = '.asc$',
                         full.names = TRUE )
# files.full # check also 
files.full[1:2]

## [1] "C:/Users/kramer/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_01.asc"
## [2] "C:/Users/kramer/PopDynIZW Dropbox/Steph Kramer/_GitHub/_DataBorneo/BaseMaps/bio_asc_02.asc"

# names only
files.rel <- list.files(path = maps_wd, pattern = '.asc$',
                        full.names = FALSE )
files.rel[1:2]

## [1] "bio_asc_01.asc" "bio_asc_02.asc"

Working with stacks

The advantage is, that you do not need to apply a command to each raster map separately, but can do it ‘all in one’. E.g., We can set the crs of each single raster in just one line. In the following, we stack four maps in a ‘geodatabase’ called ‘predictors’. Check the .doc for a description of the layers. In the following, we select 4 spatial layers:

predictors <- stack(x = files.full[c(9, 12, 22, 24)])
projection(x = predictors) # crs(predictors)

## [1] NA

projection(x = predictors) <- CRS(projargs = "+proj=longlat +datum=WGS84")

A raster stack contains the single raster layers in a list:

predictors@layers # this is a list! Address single layers with [[ ]]

## [[1]]
## class      : RasterLayer 
## dimensions : 1425, 1423, 2027775  (nrow, ncol, ncell)
## resolution : 0.008333334, 0.008333334  (x, y)
## extent     : 108.3341, 120.1925, -4.374013, 7.500988  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : bio_asc_09.asc 
## names      : bio_asc_09 
## 
## 
## [[2]]
## class      : RasterLayer 
## dimensions : 1425, 1423, 2027775  (nrow, ncol, ncell)
## resolution : 0.008333334, 0.008333334  (x, y)
## extent     : 108.3341, 120.1925, -4.374013, 7.500988  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : bio_asc_12.asc 
## names      : bio_asc_12 
## 
## 
## [[3]]
## class      : RasterLayer 
## dimensions : 1425, 1423, 2027775  (nrow, ncol, ncell)
## resolution : 0.008333334, 0.008333334  (x, y)
## extent     : 108.3341, 120.1925, -4.374013, 7.500988  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : bio_asc_24.asc 
## names      : bio_asc_24 
## values     : -2147483648, 2147483647  (min, max)
## 
## 
## [[4]]
## class      : RasterLayer 
## dimensions : 1425, 1423, 2027775  (nrow, ncol, ncell)
## resolution : 0.008333334, 0.008333334  (x, y)
## extent     : 108.3341, 120.1925, -4.374013, 7.500988  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : bio_asc_42.asc 
## names      : bio_asc_42 
## values     : -2147483648, 2147483647  (min, max)

Plotting just one layer:

plot(predictors@layers[[4]])

Plotting stacks

With base plot

plot(x = predictors)

With {tamp}: This is standardizing the legend to the min and max of all layers, which is problematic if maps are not in the same units. That’s why you cannot see the information in two of the maps:

tmap_mode(mode = "plot")
tm_shape(shp = predictors) + tm_raster()

Extract values from stacks

Violin plots to visualize summary statistics:

bwplot(x = predictors[[c(1, 2)]]) # bwplot(predictors)

# n.b. double [[ ]] because stack is a list of spatial rasters

Again, these do not look nice because of the different axis scaling.

For advanced users: More beautiful violinplots can be found in the {YaRrr} package (called pirate plots: https://bookdown.org/ndphillips/YaRrr/pirateplot.html) or simply in ggplot2. For this, we need to transform the RasterLayer data into a data.frame first:

raster_df <- na.omit(data.frame(values(predictors[[c(1,3)]])))
raster_names <- names(raster_df)
raster_ct    <- dim(raster_df)[1]
df2 <- data.frame(val = c(raster_df[,names(raster_df)[1]], 
                          raster_df[,names(raster_df)[2]]))
df2$grp <- rep(raster_names , each = raster_ct)
head(df2)

##   val        grp
## 1 274 bio_asc_09
## 2 274 bio_asc_09
## 3 275 bio_asc_09
## 4 274 bio_asc_09
## 5 274 bio_asc_09
## 6 276 bio_asc_09

#take a random subsample of the data to not crash your PC when plotting:
a <- sample(x=nrow(df2), size=1000, replace = FALSE) 
df3 <- df2[a,]

p <- ggplot(data = df3, aes(x = grp, y = val))
p + geom_violin(scale = "width",fill = "grey80", colour = "#3366FF") + 
    geom_boxplot(width = 0.2, outlier.colour = "transparent") + ## remove outliers 
    geom_jitter(height = 0, width = 0.05, alpha = 0.2, size = 1.5, colour = 'blue')

Since the units are so different, it makes more sense in this case to plot them separately, using a facet_wrap() and scales = free.

p + geom_violin(scale = "width",fill = "grey80", colour = "#3366FF") + 
    geom_boxplot(width = 0.2, outlier.colour = "transparent") + ## remove outliers 
    geom_jitter(height = 0, width = 0.05, alpha = 0.2, size = 1.5, colour = "#3366FF") +
    facet_wrap(vars(grp), scales = "free") +
    scale_x_discrete(guide = "none")  + ## remove axis ticks and labels on x
    labs(x = NULL, y = "Value") +
    theme_minimal() ## set custom plot style

# save the last plot
ggsave(paste0(output_wd, "/savedggplot.pdf"), width = 5, height = 5)
# or use here:
# ggsave(here("output", "savedggplot.pdf"), width = 5, height = 5)

…or plot them separately and combine them using the {patchwork} package:

df3_bio9 <- subset(df3, df3$grp == 'bio_asc_09')
df3_bio24 <- subset(df3, df3$grp == 'bio_asc_24')

p1 <- ggplot(data = df3_bio9, aes(x = grp, y = val)) +
      geom_violin(scale = "width",fill = "grey80", colour = "#3366FF") + 
      geom_boxplot(width = 0.2, outlier.colour = "transparent") + ## remove outliers 
      geom_jitter(height = 0, width = 0.05, alpha = 0.2, size = 1.5, colour = "#3366FF") +
      labs(x = NULL, y = "Value") +
      theme_minimal() ## set custom plot style

p2 <- ggplot(data = df3_bio24, aes(x = grp, y = val)) +
      geom_violin(scale = "width",fill = "grey80", colour = "#3366FF") + 
      geom_boxplot(width = 0.2, outlier.colour = "transparent") + ## remove outliers 
      geom_jitter(height = 0, width = 0.05, alpha = 0.2, size = 1.5, colour = "#3366FF") +
      labs(x = NULL, y = "") +
      theme_minimal() ## set custom plot style

p1 + p2

check here for how {ggplot2} works if you want to make really nice plots: * https://github.com/rstudio/cheatsheets/blob/master/data-visualization-2.1.pdf

or get inspiration here:
https://www.cedricscherer.com/2019/08/05/a-ggplot2-tutorial-for-beautiful-plotting-in-r/
https://www.r-graph-gallery.com/

Back to spatial R.

Retrieve metrics/ statistics from the rasters:

#- Extract information from all rasters in one command
# cellStats(predictors, 'mean')
round(x = cellStats(x = predictors, stat = 'mean'), digits = 2)

## bio_asc_09 bio_asc_12 bio_asc_24 bio_asc_42 
##     255.74    3013.26     251.75      10.65

Data export

Save the raster (not the plot…): There are many exchange formats for rasters. The best would be GeoTiff, which also saves the projection and is smaller than ascii. However, for modelling e.g. in MaxEnt, the raster maps are needed in ascii-format.

writeRaster(x = hillsh, 
            filename = paste0(output_wd,"/hillshade.asc"), 
            # or use here: here("output", "hillshade.asc")
            format = "ascii",
            overwrite = TRUE, 
            NAflag = -9999)

Manipulating rasters

Do any kind of raster algebra:

new_ras <- ras_bio_asc_01 + ras_bio_asc_24 + 100

For changing cell size use aggregate() or resample():

# collapse 20*20 cells into 1 using function 'mean':
ras_bio_asc_01_agg <- aggregate(x = ras_bio_asc_01, fact = 20, fun = mean)

We finally plot the two new rasters next to each other:

par(mfrow = c(1,2)) # recap from course 1: plots in 1 row and 2 columns
plot(x = new_ras, col = rev(rainbow(5)))
plot(ras_bio_asc_01_agg)

par(mfrow = c(1,1)) # remember to set it back to the default

Change and query raster values

# convert ras_bio_asc_01 to degree Celsius units (divide by 10)
range(values(ras_bio_asc_01), na.rm = TRUE)

## [1]  63 278

mean.t.c <- ras_bio_asc_01 / 10
range(values(mean.t.c), na.rm = TRUE)

## [1]  6.3 27.8

# find areas with mean annual temp >= 25 deg C
mean.t.c.25 <- mean.t.c >= 25

Have a look at the result using a binary plot:

plot(x = mean.t.c.25)

Accessing cell values

Remember the function extract(). In the following, we extract from the small raster the first three values (= cols 1-3; longitude or x-axis) at line 5 (= row 5; latitude or y-axis).

cells <- cellFromRowCol(object = ras_bio_asc_01_cr, row = 5, col = 1:3)
cells    #gives cell ID number, the index/ rownumber!!!

## [1] 121 122 123

extract(x = ras_bio_asc_01_cr, y = cells) # gives cell values!

## [1] 266 265 266

plot(x = ras_bio_asc_01_cr)
# to plot the points on top, insert the 'cells' index into 
# the dataframe of the coordinates of the RasterLayer object ras_bio_asc_01_cr
points(x = coordinates(obj = ras_bio_asc_01_cr)[cells,], col = "blue")

Extract rastervalues from all rasters at once (raster stack):

This is very important if you want to create your master table, i.e. based on the xy-coordinates of your species sightings (i.e. GPS point locations), you could extract the important environmental predictors at that location.

We do that for an area in the center of Borneo. For this, we first get the row and column indices of the center of the large Borneo map, e.g. ras_bio_asc_01, and 5 x 5 cells in addition, i.e. a square with an area of 5 km x 5 km:

center_x = floor(nrow(predictors)/2) # learn about the functions round(), 
center_y = floor(ncol(predictors)/2) #                 ceiling(), floor()
center_x; center_y

## [1] 712

## [1] 711

stack_cells <- cellFromRowCol(
  object = ras_bio_asc_01, 
  row = center_x:(center_x + 5), 
  col = center_y:(center_y + 5)
)
head(stack_cells)

## [1] 1012464 1013888 1015312 1016736 1018160 1019584

Why are these index numbers so large? Because the counting starts at the upper left corner of the map and continues consecutively, i.e. the data (ras_bio_asc_01@data@values) are one huge vector with length nrow * ncol!

# finally, extract them:
pred_dat <- extract(predictors, stack_cells) 
head(pred_dat)

##      bio_asc_09 bio_asc_12 bio_asc_24 bio_asc_42
## [1,]        254       4276        303          1
## [2,]        249       4206        352          1
## [3,]        243       4129        439          1
## [4,]        232       3999        702          2
## [5,]        232       4003        588          2
## [6,]        218       3810        881          2

mastertable <- data.frame(stack_cells, pred_dat)
mastertable

##   stack_cells bio_asc_09 bio_asc_12 bio_asc_24 bio_asc_42
## 1     1012464        254       4276        303          1
## 2     1013888        249       4206        352          1
## 3     1015312        243       4129        439          1
## 4     1016736        232       3999        702          2
## 5     1018160        232       4003        588          2
## 6     1019584        218       3810        881          2

Change raster values:

xy = xyFromCell(object = ras_bio_asc_01_cr, cell = cells) # coordinates(ras_bio_asc_01_cr)[cells,]
extract(x = ras_bio_asc_01_cr, y = xy)          # extract(ras_bio_asc_01_cr, cells)

## [1] 266 265 266

# Change values, e.g. for adding forest or creating a corridor
# take care! irreversible! better work on a copy!
copy_ras <- ras_bio_asc_01_cr
copy_ras[cells] <- 250
plot(x = copy_ras)

Compute distance to points

# calculate distance
my_dist <- distanceFromPoints(object = ras_bio_asc_01_cr, xy = xy)
plot(x = my_dist) #units?

# nicer plot adding the points from which a distance should be computed.
xy_sf <- st_as_sf(x = data.frame(xy), 
                  coords = c("x","y"),
                  crs = 4326,
                  sf_column_name = "geometry" ) 

tmap_mode(mode = "plot")
tm_shape(shp = my_dist) + 
    tm_raster(n = 100, palette = rev(grDevices::terrain.colors(100)),
              legend.show = FALSE) +
  tm_shape(shp = xy_sf) + 
  tm_dots(size = 1)

Calculate distance from many points:

cells1 <- c(cells, 250, 360) # add two more points
xy = xyFromCell(ras_bio_asc_01_cr, cells1) # =coordinates(Bor_mat_cr)[cells,]
my_dist <- distanceFromPoints(ras_bio_asc_01_cr, xy)
plot(my_dist) #units?

Vector data / shapefiles

Polygons and lines

Import shapefiles

Currently, there are two packages sp and sf (standing for simple features), both with still important functionality. However, sf is much easier to use and handle. The suggestion currently is: learn both!

In the following, the most important commands are provided for the sf-package, and if necessary, also for the sp-package:

# Border of countries and provinces of Borneo
# - only loading columns 1:3, 5, 7, 17, 18 of attribute table:
Borneo_shp <- st_read(dsn = maps_wd,
                      layer = "borneo_admin",
                      stringsAsFactors = FALSE)[, c(1:3, 5, 7, 17, 18)]

## Reading layer `borneo_admin' from data source 
##   `C:\Users\kramer\PopDynIZW Dropbox\Steph Kramer\_GitHub\_DataBorneo\BaseMaps' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 158 features and 18 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 108.5986 ymin: -4.943054 xmax: 119.2697 ymax: 7.380556
## Geodetic CRS:  WGS 84

# Protected areas (National Parks, Nature Reserves, Forest Reserves)
PA_shp     <-  st_read(dsn = maps_wd,
                       layer = "Bor_PA",
                       stringsAsFactors = FALSE)[, c(1:4)]

## Reading layer `Bor_PA' from data source 
##   `C:\Users\kramer\PopDynIZW Dropbox\Steph Kramer\_GitHub\_DataBorneo\BaseMaps' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 255 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 108.5969 ymin: -4.18356 xmax: 119.6176 ymax: 9.911229
## Geodetic CRS:  WGS 84

# fix problematic polygons
PA_shp <- st_make_valid(PA_shp)

# Main rivers
River_shp  <- st_read(dsn = maps_wd,
                      layer = "sn_100000",
                      stringsAsFactors = FALSE)

## Reading layer `sn_100000' from data source 
##   `C:\Users\kramer\PopDynIZW Dropbox\Steph Kramer\_GitHub\_DataBorneo\BaseMaps' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 396 features and 4 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 108.9454 ymin: -3.787917 xmax: 118.7879 ymax: 6.464583
## Geodetic CRS:  WGS 84

## the 'old' sp package would work with that command:
Admin_shp <- readOGR(dsn = maps_wd, layer = "borneo_admin",
                     stringsAsFactors = FALSE)[,c(1:3,5,7,17,18)]

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\kramer\PopDynIZW Dropbox\Steph Kramer\_GitHub\_DataBorneo\BaseMaps", layer: "borneo_admin"
## with 158 features
## It has 18 fields

## transformations
Admin_sf <- as(Admin_shp, "sf") ## from sp object to sf object 
Admin_sp <- as(Admin_sf, "Spatial") ## from sf to sp object

Working with vector data

Accessing simple feature objects — easy handling, similar to data.frames:

#Please note the similarity to accessing info from rasters.
str(object = extent(PA_shp))

## Formal class 'Extent' [package "raster"] with 4 slots
##   ..@ xmin: Named num 109
##   .. ..- attr(*, "names")= chr "xmin"
##   ..@ xmax: Named num 120
##   .. ..- attr(*, "names")= chr "xmax"
##   ..@ ymin: Named num -4.18
##   .. ..- attr(*, "names")= chr "ymin"
##   ..@ ymax: Named num 9.91
##   .. ..- attr(*, "names")= chr "ymax"

extent(PA_shp)@xmin

##     xmin 
## 108.5969

names(x = PA_shp) #returns column names of a.t.

## [1] "SITE_ID"  "NAME_ENG" "COUNTRY"  "SUB_NAT"  "geometry"

Accessing SpatialPolygonsDataFrames:

# Much more complex than sf-objects.
class(Admin_shp)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

str(object = extent(Admin_shp))

## Formal class 'Extent' [package "raster"] with 4 slots
##   ..@ xmin: num 109
##   ..@ xmax: num 119
##   ..@ ymin: num -4.94
##   ..@ ymax: num 7.38

extent(Admin_shp)@xmin

## [1] 108.5986

names(x = Admin_shp) #returns column names of a.t.

## [1] "ID_0"       "ISO"        "NAME_0"     "NAME_1"     "NAME_2"    
## [6] "Shape_Leng" "Shape_Area"

Summarising spatial information for each column of attribute table [a.t.]:

summary(object = PA_shp)

##     SITE_ID         NAME_ENG           COUNTRY            SUB_NAT         
##  Min.   :   783   Length:255         Length:255         Length:255        
##  1st Qu.:  3844   Class :character   Class :character   Class :character  
##  Median :  9842   Mode  :character   Mode  :character   Mode  :character  
##  Mean   : 45309                                                           
##  3rd Qu.: 19582                                                           
##  Max.   :901238                                                           
##           geometry  
##  MULTIPOLYGON :255  
##  epsg:4326    :  0  
##  +proj=long...:  0  
##                     
##                     
##

str(object = PA_shp)

## Classes 'sf' and 'data.frame':   255 obs. of  5 variables:
##  $ SITE_ID : int  785 787 791 795 1300 783 784 786 788 789 ...
##  $ NAME_ENG: chr  "Kinabalu" "Mulu" "Niah" "Tawau Hill Park" ...
##  $ COUNTRY : chr  "Malaysia" "Malaysia" "Malaysia" "Malaysia" ...
##  $ SUB_NAT : chr  "Sabah" "Sarawak" "Sarawak" "Sabah" ...
##  $ geometry:sfc_MULTIPOLYGON of length 255; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:322, 1:2] 117 117 117 117 117 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA
##   ..- attr(*, "names")= chr [1:4] "SITE_ID" "NAME_ENG" "COUNTRY" "SUB_NAT"

PA_shp

## Simple feature collection with 255 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 108.5969 ymin: -4.18356 xmax: 119.6176 ymax: 9.911229
## Geodetic CRS:  WGS 84
## First 10 features:
##    SITE_ID        NAME_ENG  COUNTRY SUB_NAT                       geometry
## 1      785        Kinabalu Malaysia   Sabah MULTIPOLYGON (((116.5979 6....
## 2      787            Mulu Malaysia Sarawak MULTIPOLYGON (((114.9033 4....
## 3      791            Niah Malaysia Sarawak MULTIPOLYGON (((113.7986 3....
## 4      795 Tawau Hill Park Malaysia   Sabah MULTIPOLYGON (((117.8497 4....
## 5     1300  Lanjak-Entimau Malaysia Sarawak MULTIPOLYGON (((112.2054 1....
## 6      783        Semporna Malaysia   Sabah MULTIPOLYGON (((118.7796 4....
## 7      784        Samunsam Malaysia Sarawak MULTIPOLYGON (((109.6615 1....
## 8      786       Similajau Malaysia Sarawak MULTIPOLYGON (((113.1826 3....
## 9      788           Klias Malaysia   Sabah MULTIPOLYGON (((115.636 5.3...
## 10     789      Pulau Tiga Malaysia   Sabah MULTIPOLYGON (((115.6796 5....

attributes(x = PA_shp)

## $names
## [1] "SITE_ID"  "NAME_ENG" "COUNTRY"  "SUB_NAT"  "geometry"
## 
## $row.names
##   [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
##  [19]  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36
##  [37]  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54
##  [55]  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72
##  [73]  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90
##  [91]  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108
## [109] 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
## [127] 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
## [145] 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162
## [163] 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## [181] 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
## [199] 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
## [217] 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
## [235] 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252
## [253] 253 254 255
## 
## $sf_column
## [1] "geometry"
## 
## $agr
##  SITE_ID NAME_ENG  COUNTRY  SUB_NAT 
##     <NA>     <NA>     <NA>     <NA> 
## Levels: constant aggregate identity
## 
## $class
## [1] "sf"         "data.frame"

Have a look at the content:

head(x = PA_shp)

## Simple feature collection with 6 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 111.8754 ymin: 1.32105 xmax: 118.7899 ymax: 6.4952
## Geodetic CRS:  WGS 84
##   SITE_ID        NAME_ENG  COUNTRY SUB_NAT                       geometry
## 1     785        Kinabalu Malaysia   Sabah MULTIPOLYGON (((116.5979 6....
## 2     787            Mulu Malaysia Sarawak MULTIPOLYGON (((114.9033 4....
## 3     791            Niah Malaysia Sarawak MULTIPOLYGON (((113.7986 3....
## 4     795 Tawau Hill Park Malaysia   Sabah MULTIPOLYGON (((117.8497 4....
## 5    1300  Lanjak-Entimau Malaysia Sarawak MULTIPOLYGON (((112.2054 1....
## 6     783        Semporna Malaysia   Sabah MULTIPOLYGON (((118.7796 4....

tail(x = PA_shp)

## Simple feature collection with 6 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 111.4107 ymin: -3.035619 xmax: 118.2907 ymax: 6.4952
## Geodetic CRS:  WGS 84
##     SITE_ID                NAME_ENG           COUNTRY SUB_NAT
## 250  317300        Maganting Island          Malaysia   Sabah
## 251  317301                 Bod Tai          Malaysia   Sabah
## 252  342562 Tanjung Penghujan NR/RP         Indonesia    <NA>
## 253  901236          Tasek Merimbun Brunei Darussalam    <NA>
## 254  901237  Kinabalu National Park          Malaysia    <NA>
## 255  901238      Mulu National Park          Malaysia    <NA>
##                           geometry
## 250 MULTIPOLYGON (((118.2839 4....
## 251 MULTIPOLYGON (((118.2182 5....
## 252 MULTIPOLYGON (((111.4377 -2...
## 253 MULTIPOLYGON (((114.7238 4....
## 254 MULTIPOLYGON (((116.5979 6....
## 255 MULTIPOLYGON (((114.9033 4....

Retrieving information of a.t. — recap working with data.frames from Day1_R-Intro course.

PA_shp[1, ] #returns first entry (row) of all 4 columns

## Simple feature collection with 1 feature and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 116.4638 ymin: 5.99952 xmax: 116.7631 ymax: 6.4952
## Geodetic CRS:  WGS 84
##   SITE_ID NAME_ENG  COUNTRY SUB_NAT                       geometry
## 1     785 Kinabalu Malaysia   Sabah MULTIPOLYGON (((116.5979 6....

PA_shp[, 2] # returns summary of col only!

## Simple feature collection with 255 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 108.5969 ymin: -4.18356 xmax: 119.6176 ymax: 9.911229
## Geodetic CRS:  WGS 84
## First 10 features:
##           NAME_ENG                       geometry
## 1         Kinabalu MULTIPOLYGON (((116.5979 6....
## 2             Mulu MULTIPOLYGON (((114.9033 4....
## 3             Niah MULTIPOLYGON (((113.7986 3....
## 4  Tawau Hill Park MULTIPOLYGON (((117.8497 4....
## 5   Lanjak-Entimau MULTIPOLYGON (((112.2054 1....
## 6         Semporna MULTIPOLYGON (((118.7796 4....
## 7         Samunsam MULTIPOLYGON (((109.6615 1....
## 8        Similajau MULTIPOLYGON (((113.1826 3....
## 9            Klias MULTIPOLYGON (((115.636 5.3...
## 10      Pulau Tiga MULTIPOLYGON (((115.6796 5....

# PA_shp$NAME_ENG #returns a long list  
head(x = PA_shp$NAME_ENG) #using fct head() to only show the first 6 entries

## [1] "Kinabalu"        "Mulu"            "Niah"            "Tawau Hill Park"
## [5] "Lanjak-Entimau"  "Semporna"

PA_shp[1, 1] # = PA_shp$SITE_ID[1]

## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 116.4638 ymin: 5.99952 xmax: 116.7631 ymax: 6.4952
## Geodetic CRS:  WGS 84
##   SITE_ID                       geometry
## 1     785 MULTIPOLYGON (((116.5979 6....

PA_shp[2, 3] # = PA_shp$COUNTRY[2]

## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 114.7777 ymin: 3.9186 xmax: 114.9945 ymax: 4.28586
## Geodetic CRS:  WGS 84
##    COUNTRY                       geometry
## 2 Malaysia MULTIPOLYGON (((114.9033 4....

The following returns the row indices of the data frame or a.t., respectively:

which(x = PA_shp$COUNTRY == 'Malaysia')

##   [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  26  27  28  29  30
##  [19]  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48
##  [37]  49  50  51  52  53  54  55  56  63  64  65  66  67  68  69  70  71  72
##  [55]  73  74  75  76  77  80  81 101 102 103 104 113 114 115 116 117 118 119
##  [73] 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
##  [91] 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155
## [109] 156 157 158 159 161 162 163 164 165 166 167 169 170 171 172 185 186 187
## [127] 190 191 192 193 194 195 196 221 222 223 224 225 226 227 228 230 233 237
## [145] 238 239 240 241 249 250 251 254 255

Visualising vector data

You can use many packages, e.g. {ggplot2}, for {sf} objects: https://www.r-spatial.org/r/2018/10/25/ggplot2-sf.html.

However, it does NOT work for {sp} objects (SpatialPolygonsDataFrame).

ggplot(data = PA_shp) +
  geom_sf(fill = "chartreuse3") + 
  labs(
    x = "Longitude", y =  "Latitude",
    title = "Protected areas"
  ) +
  theme_minimal() ## set custom plot style

or {tamp}:

tmap_mode(mode = "plot")
tm_shape(shp = PA_shp[, 1]) + 
  tm_polygons(col = "SITE_ID", palette = grDevices::terrain.colors(5),
              style = "quantile") # color scheme by quantiles

Interactive mode:

tmap_mode(mode = "view")
tm_shape(shp = PA_shp[, 1]) + 
  tm_polygons(col = "SITE_ID", palette = grDevices::terrain.colors(5),
              style = "quantile") # color scheme  by quantiles

It is also possible to use a simple plot() function for {sf} objects. However, mind that this object has 7 columns, so each column holds information that will be plotted, i.e. you will plot 7 single plots. So do a little transformation: either plot just one column, or define it as spatial:

# plot(Borneo_shp) not run
plot(Borneo_shp[,1]) # select column

plot(as(Borneo_shp, "Spatial")) #plot geometry without information of a col

And for the {sp} package:

Here, the basic plot function works:

plot(Admin_sp)

There is also a special function for plotting {sp} objects, however, also here, you need to select which information you want to plot.

# spplot(Admin_sp) #Don't do that! Each column of the a.t. will be plotted
head(Admin_sp)

##   ID_0 ISO   NAME_0 NAME_1     NAME_2 Shape_Leng Shape_Area
## 0  158 MYS Malaysia  Sabah Lahad Datu  6.7475367 0.62106012
## 1  158 MYS Malaysia  Sabah      Papar  1.6727672 0.10065591
## 2  158 MYS Malaysia  Sabah  Penampang  0.9548857 0.03951545
## 3  158 MYS Malaysia  Sabah Pensiangan  4.0257859 0.49725547
## 4  158 MYS Malaysia  Sabah      Pitas  2.9223456 0.11955352
## 5  158 MYS Malaysia  Sabah      Ranau  2.4815704 0.24027591

spplot(Admin_sp[6]) # only plot one of the geometries, coloured e.g. by shape_length

Plot a detail on the map — plot of a single selected protected area:

tmap_mode(mode = "plot")
tm_shape(shp = Borneo_shp) + 
    tm_polygons(border.col = "deepskyblue4") +
  tm_shape(PA_shp[,1]) + 
    tm_polygons(border.col = "black") +
  tm_shape(PA_shp[1, ]) + 
    tm_polygons(border.col = "red")

Small exercises

Simplify your life! Work with subsets, select Malaysia, save and plot as shapefile:

Mal_PA_shp <- subset(PA_shp, PA_shp$COUNTRY == 'Malaysia')
st_write(obj = Mal_PA_shp, 
         dsn = output_wd, 
         layer = 'test',
         driver = 'ESRI Shapefile', 
         delete_layer = TRUE)

## Deleting layer `test' using driver `ESRI Shapefile'
## Writing layer `test' to data source 
##   `C:/Users/kramer/PopDynIZW Dropbox/Steph Kramer/_GitHub/output' using driver `ESRI Shapefile'
## Writing 153 features with 4 fields and geometry type Multi Polygon.

tmap_mode(mode = "plot")
tm_shape(shp = Borneo_shp) + 
    tm_polygons(border.col = "blue") +
  tm_shape(shp = Mal_PA_shp) + 
    tm_polygons(col = "green")

Extract elevation per PA and save average to a.t. as a new column

fewPA <- Mal_PA_shp[c(1:5), 1]
tmp <- extract(x = ras_bio_asc_24, y = fewPA)   # returns a list — each element contains
                                       # the elevation raster cells (ras_bio_asc_24)
str(tmp)

## List of 5
##  $ : num [1:884] 286 491 508 500 775 986 645 582 335 675 ...
##  $ : num [1:652] 75 68 63 252 121 63 66 53 49 249 ...
##  $ : num [1:40] 6 18 29 24 30 23 19 18 29 17 ...
##  $ : num [1:326] 182 135 168 269 138 214 327 444 107 138 ...
##  $ : num [1:1779] 120 112 148 176 202 158 111 176 153 122 ...

#lapply(tmp, FUN = summary)

# Append mean elevation:
fewPA$mean_elevation <- round(x = unlist(lapply(tmp, FUN = mean, na.rm = TRUE)))
fewPA

## Simple feature collection with 5 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 111.8754 ymin: 1.32105 xmax: 118.0747 ymax: 6.4952
## Geodetic CRS:  WGS 84
##   SITE_ID                       geometry mean_elevation
## 1     785 MULTIPOLYGON (((116.5979 6....           1348
## 2     787 MULTIPOLYGON (((114.9033 4....            555
## 3     791 MULTIPOLYGON (((113.7986 3....            108
## 4     795 MULTIPOLYGON (((117.8497 4....            512
## 5    1300 MULTIPOLYGON (((112.2054 1....            404

Geospatial calculation

Using the ‘old’ {sp} package and the {rgeos} package, commands for calculating properties of spatial objects were: gArea(), gLength(), gDistance(). Other often used GIS functions relate to basic spatial operations, such as gBuffer(), gCentroid(), gDifference(), gUnion(), and gIntersection(), gUnionCascaded() (= dissolve function) and gSimplify(). Spatial queries were gIntersects(), gContains(), gIsValid(). http://www.nickeubank.com/wp-content/uploads/2015/10/RGIS2_MergingSpatialData_part2_GeometricManipulations.html.

The new {sf} package has the same functionality, with functions starting with st_, e.g. st_buffer(), st_intersection() etc.
https://cran.r-project.org/web/packages/sf/vignettes/sf1.html
https://geocompr.robinlovelace.net/spatial-operations.html

Spatial query

Returns the area in square meters:

# st_area(x = Borneo_shp) #returns long vector
head(st_area(x = Borneo_shp))

## Units: [m^2]
## [1] 7651059265 1238543644  486042067 6128271711 1468016456 2955188183

Spatial transformation (projection)

Borneo_shp_moll <-  st_transform(Borneo_shp, c("+proj=moll +datum=WGS84"))
class(Borneo_shp_moll) #sf object, data.frame!

## [1] "sf"         "data.frame"

tmap_mode(mode = "plot")
tm_shape(shp = Borneo_shp_moll) + tm_polygons(border.col = "blue") #do you see the difference?

head(st_area(x = Borneo_shp_moll)) # units?

## Units: [m^2]
## [1] 7666715202 1241315346  487129643 6141990607 1471302821 2961802507

ggplot(Borneo_shp_moll) +
  geom_sf(fill = "chartreuse3") +
  theme_minimal() ## set custom plot style

Example: Area of Malaysia in Borneo.

Mal_Borneo_shp <- subset(Borneo_shp_moll, Borneo_shp_moll$NAME_0 == 'Malaysia')
head(st_area(x = Mal_Borneo_shp) / 1000000) # or / 1e6

## Units: [m^2]
## [1] 7666.7152 1241.3153  487.1296 6141.9906 1471.3028 2961.8025

# or use set_units to change the units from m^2 to km^2
head(set_units(x = st_area(x = Mal_Borneo_shp), value = km^2))

## Units: [km^2]
## [1] 7666.7152 1241.3153  487.1296 6141.9906 1471.3028 2961.8025

Example: Calculate area of a polygon and add the area of a polygon as new column to a.t.

st_area(x = Borneo_shp_moll[3, ]) # for a single polygon

## 487129643 [m^2]

head(st_area(x = Borneo_shp_moll, byid = TRUE)) #for all polygons

## Units: [m^2]
## [1] 7666715202 1241315346  487129643 6141990607 1471302821 2961802507

area_km2 <- set_units(x = st_area(x = Borneo_shp_moll, byid = TRUE), value = km^2)
Borneo_shp_moll = data.frame(Borneo_shp_moll, area_km2)
head(x = Borneo_shp_moll)

##   ID_0 ISO   NAME_0 NAME_1     NAME_2 Shape_Leng Shape_Area
## 1  158 MYS Malaysia  Sabah Lahad Datu  6.7475367 0.62106012
## 2  158 MYS Malaysia  Sabah      Papar  1.6727672 0.10065591
## 3  158 MYS Malaysia  Sabah  Penampang  0.9548857 0.03951545
## 4  158 MYS Malaysia  Sabah Pensiangan  4.0257859 0.49725547
## 5  158 MYS Malaysia  Sabah      Pitas  2.9223456 0.11955352
## 6  158 MYS Malaysia  Sabah      Ranau  2.4815704 0.24027591
##                         geometry         area_km2
## 1 MULTIPOLYGON (((11824461 58... 7666.7152 [km^2]
## 2 MULTIPOLYGON (((11590069 72... 1241.3153 [km^2]
## 3 MULTIPOLYGON (((11592912 72...  487.1296 [km^2]
## 4 MULTIPOLYGON (((11651495 62... 6141.9906 [km^2]
## 5 MULTIPOLYGON (((11703665 83... 1471.3028 [km^2]
## 6 MULTIPOLYGON (((11675075 77... 2961.8025 [km^2]

The following for {sp} objects:

#gArea(Borneo_shp_moll, byid=TRUE)/1e6 # does not work, as Borneo_shp_moll is sf-object!
head(gArea(Admin_sp, byid=TRUE)/1e6) # what does the warning message mean? -> check crs!

##            0            1            2            3            4            5 
## 6.210601e-07 1.006559e-07 3.951545e-08 4.972555e-07 1.195535e-07 2.402759e-07

Admin_sp_moll <- spTransform(Admin_sp, c("+proj=moll +datum=WGS84"))
head(gArea(Admin_sp_moll, byid=TRUE)/1e6) # now it works!

##         0         1         2         3         4         5 
## 7666.7152 1241.3153  487.1296 6141.9906 1471.3028 2961.8025

Points

Import spatial points

Usually, shapefiles come with a .proj file defining the projection, but not always. Please always check whether the projection is set:

pt_shp <- st_read(dsn = paste0(recs_wd, "/FCsim.shp"))

## Reading layer `FCsim' from data source 
##   `C:\Users\kramer\PopDynIZW Dropbox\Steph Kramer\_GitHub\_DataBorneo\BaseRecords\FCsim.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 55 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 109.3758 ymin: -2.744159 xmax: 118.4117 ymax: 5.824283
## CRS:           NA

pt_shp  # crs is missing!

## Simple feature collection with 55 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 109.3758 ymin: -2.744159 xmax: 118.4117 ymax: 5.824283
## CRS:           NA
## First 10 features:
##    Id  POINT_X  POINT_Y                  geometry
## 1   0 118.4117 5.221667 POINT (118.4117 5.221667)
## 2   0 118.2645 5.005603 POINT (118.2645 5.005603)
## 3   0 117.7992 5.029167 POINT (117.7992 5.029167)
## 4   0 117.6058 5.177957 POINT (117.6058 5.177957)
## 5   0 117.5566 5.381088 POINT (117.5566 5.381088)
## 6   0 117.4704 5.005603 POINT (117.4704 5.005603)
## 7   0 117.2734 5.227201 POINT (117.2734 5.227201)
## 8   0 116.7256 5.319533 POINT (116.7256 5.319533)
## 9   0 116.6467 5.162500   POINT (116.6467 5.1625)
## 10  0 117.2475 4.888333 POINT (117.2475 4.888333)

st_crs(pt_shp) = 4326 # set it with command st_crs
pt_shp

## Simple feature collection with 55 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 109.3758 ymin: -2.744159 xmax: 118.4117 ymax: 5.824283
## Geodetic CRS:  WGS 84
## First 10 features:
##    Id  POINT_X  POINT_Y                  geometry
## 1   0 118.4117 5.221667 POINT (118.4117 5.221667)
## 2   0 118.2645 5.005603 POINT (118.2645 5.005603)
## 3   0 117.7992 5.029167 POINT (117.7992 5.029167)
## 4   0 117.6058 5.177957 POINT (117.6058 5.177957)
## 5   0 117.5566 5.381088 POINT (117.5566 5.381088)
## 6   0 117.4704 5.005603 POINT (117.4704 5.005603)
## 7   0 117.2734 5.227201 POINT (117.2734 5.227201)
## 8   0 116.7256 5.319533 POINT (116.7256 5.319533)
## 9   0 116.6467 5.162500   POINT (116.6467 5.1625)
## 10  0 117.2475 4.888333 POINT (117.2475 4.888333)

Spatial points are defined by two columns with x and y coordinates, can be imported as normal data.frames and then converted to spatial objects:

pt_file <- paste0(recs_wd, "/MyNewSpecies.csv")
df_recs <- read.table(file = pt_file, header = TRUE, sep = ',')
class(x = df_recs)

## [1] "data.frame"

head(x = df_recs)

##   species     long         lat
## 1 mysimsp 109.4716 -1.11984615
## 2 mysimsp 109.3216 -0.02817942
## 3 mysimsp 117.0383  3.60515411
## 4 mysimsp 109.1550  0.09682059
## 5 mysimsp 111.1550  1.78848735
## 6 mysimsp 115.4550  4.93848752

See below of how to convert a data.frame into a spatial {sf} object.

Plot the points

Simple data frames with location information, as long as the coordinates are the same, can also be plotted in space. We first plot the polygon:

plot(x = as(Borneo_shp, "Spatial"), col = 'grey', border = 'white') # polygon
points(x = df_recs$long, df_recs$lat, cex = 0.5, pch = 15) #simple d.f.! 
plot(x = as(pt_shp, "Spatial"), col = 'blue', add = TRUE)

{tmap} can only handle spatial objects, so the dataframe df_recs needs to be converted into an {sf} object first. You have to define which columns contain the coordinates:

recs_sf <- st_as_sf(x = data.frame(df_recs),
                    coords = c("long", "lat"),
                    crs = 4326,
                    sf_column_name = "geometry")

Plot it with {tamp}:

tmap_mode(mode = "plot")
tm_shape(shp = Borneo_shp) + 
    tm_polygons(col = "grey", border.col = "white") +
  tm_shape(shp = recs_sf) + 
    tm_dots(shape = 3, size = 0.5) +
  tm_shape(shp = pt_shp) + 
    tm_dots(col = "blue")

Spatial overlay — vector data

# retrieve the geometry (location) indices of PA_shp at
# the locations of sp_recs: which points are in PA_shp
# insidePA <- st_intersection(x = recs_sf, y = PA_shp)  ## doesn't run, fix polygons first
insidePA <- st_intersection(x = recs_sf, y = PA_shp)

Spatial overlay — vector and raster

# for a RASTER: extract mean ann. temp. from ras_bio_asc_01
# and add it to a.t.
mean_t <- extract(x = ras_bio_asc_01, y = recs_sf)
recs_sf$mean_t <- mean_t
mean(x = recs_sf$mean_t) # hist(sp_recs_sf$mean_t)

## [1] 272.7788

Saving vector data

st_write(obj = insidePA,
         dsn = output_wd, 
         layer = "inPA",
         driver = "ESRI Shapefile",
         delete_layer = TRUE)

## Deleting layer `inPA' using driver `ESRI Shapefile'
## Writing layer `inPA' to data source 
##   `C:/Users/kramer/PopDynIZW Dropbox/Steph Kramer/_GitHub/output' using driver `ESRI Shapefile'
## Writing 11 features with 5 fields and geometry type Point.

At the end, a beautiful plot:

tmap_mode(mode = "plot")
tm_shape(shp = hillsh) + 
  tm_raster(palette = "Greys",
            legend.show = FALSE) +
  tm_shape(shp = ras_bio_asc_24) + 
  tm_raster(palette = terrain.colors(25),
            alpha = 0.3,
            legend.show = FALSE) +
  tm_shape(shp = River_shp) + 
  tm_lines(col = "dodgerblue3") +
  tm_shape(shp = recs_sf) + 
  tm_dots(shape = 16, 
          size = 0.1) +
  tm_shape(shp = insidePA) + 
  tm_dots(col = "red",
          size = 0.5,
          shape = 3) +
  tm_shape(shp = PA_shp) + 
  tm_polygons(border.col = "white", 
              alpha = 0) +
  tm_layout(title = "END OF SESSION", 
            title.color = "darkgrey",
            title.position = c(0.01,0.9),
            title.size = 2)

Exercise 1

Revisit the data set from course 1 on wild boar observations in Berlin: data_wb_melden_en.csv. Now, we want to study the spatial patterns of the wild boar observations. We may hypothesize that wild boar observations are also related to local differences in weather. Answer the following question using spatial data sets and visualizations:

Question 1.1) Were wild boar with many piglets seen more often in warm parts of the city?
Question 1.2) Were wild boars more often observed at colder spots of the city during sunny weather?
Question 1.3) Additional question for the fast ones: Were large groups of wild boar more frequently seen in areas providing a dense tree cover?

Follow these steps to answer the questions:

1. Load the wild boar data (data_wb_melden_en.csv) and spatial data on temperature in Berlin (summer_temp_100m_3035.asc in _DataBerlin). Remember to set the correct CRS to the raster! (Hint: we need the Lambert Azimuthal Equal Area (LAEA) projection for Europe. Tip: always save the CRS to the filename ;-) ).
1. Turn the wild boar data into a simple features object. Remember to set the correct CRS (the wild boar locations were collected in WGS 84) and to transform it to the same CRS as the raster afterwards!
1. Inspect both spatial data sets by plotting the temperature raster with the wild boar locations on top using the R basic plot function. The locations should match the map (if not something went wrong when setting the CRS; check if they are the same).
1. Select one of the temperature layers and extract the values for each wild boar location.
1. Create a plot that visualizes the temperature at each wild boar observation (y axis) for each weather category (e.g. box plot or violin plot + raw data).
1. Create a new variable that holds the weather category as either “sunny” or “other”.
1. Visualize the temperature at each wild boar observation (y axis) for sunny and other weather (box and whisker plot + raw data).
1. Save the wild boar observations with the local temperature information as an .Rds file.

For the additional question:

1. Load the tree cover data set (tree_cover_density_2012_100m_3035.asc in _DataBerlin).
1. Inspect the data set by plotting the raster using the R basic plot function.
1. Extract the tree cover within a buffer of 100m around each wild boar location (hint: use the help to inspect the arguments of extract()).
1. Create a plot that shows the tree cover (y axis) based on the group size (box plot + raw data).

END

Adding the session info can be very helpful when going back to old scripts or using scripts of others:

sessionInfo()

## R version 4.1.2 (2021-11-01)
## Platform: i386-w64-mingw32/i386 (32-bit)
## Running under: Windows 10 x64 (build 19043)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_Germany.1252  LC_CTYPE=English_Germany.1252   
## [3] LC_MONETARY=English_Germany.1252 LC_NUMERIC=C                    
## [5] LC_TIME=C                       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] patchwork_1.1.1    ggplot2_3.3.5      units_0.7-2        viridis_0.6.2     
##  [5] viridisLite_0.4.0  GISTools_0.7-4     MASS_7.3-54        RColorBrewer_1.1-2
##  [9] maptools_1.1-2     dismo_1.3-5        tmap_3.3-2         stars_0.5-5       
## [13] abind_1.4-5        sf_1.0-4           rasterVis_0.51.1   lattice_0.20-45   
## [17] rgl_0.108.3        rgeos_0.5-9        rgdal_1.5-27       raster_3.5-11     
## [21] sp_1.4-6           here_1.0.1        
## 
## loaded via a namespace (and not attached):
##  [1] rprojroot_2.0.2         tools_4.1.2             bslib_0.3.1            
##  [4] utf8_1.2.2              R6_2.5.1                KernSmooth_2.23-20     
##  [7] DBI_1.1.1               colorspace_2.0-2        withr_2.4.3            
## [10] tidyselect_1.1.1        gridExtra_2.3           leaflet_2.0.4.1        
## [13] compiler_4.1.2          textshaping_0.3.6       leafem_0.1.6           
## [16] labeling_0.4.2          bookdown_0.24           sass_0.4.0             
## [19] scales_1.1.1            classInt_0.4-3          hexbin_1.28.2          
## [22] proxy_0.4-26            systemfonts_1.0.3       stringr_1.4.0          
## [25] digest_0.6.29           foreign_0.8-81          rmarkdown_2.11         
## [28] base64enc_0.1-3         dichromat_2.0-0         jpeg_0.1-9             
## [31] pkgconfig_2.0.3         htmltools_0.5.2         fastmap_1.1.0          
## [34] highr_0.9               htmlwidgets_1.5.4       rlang_0.4.12           
## [37] jquerylib_0.1.4         generics_0.1.1          farver_2.1.0           
## [40] zoo_1.8-9               jsonlite_1.7.3          crosstalk_1.2.0        
## [43] dplyr_1.0.7             magrittr_2.0.1          s2_1.0.7               
## [46] Rcpp_1.0.8              munsell_0.5.0           fansi_1.0.2            
## [49] lifecycle_1.0.1         terra_1.4-22            stringi_1.7.6          
## [52] leafsync_0.1.0          yaml_2.2.1              tmaptools_3.1-1        
## [55] grid_4.1.2              parallel_4.1.2          crayon_1.4.2           
## [58] knitr_1.36              pillar_1.6.4            codetools_0.2-18       
## [61] wk_0.5.0                XML_3.99-0.8            glue_1.6.0             
## [64] evaluate_0.14           leaflet.providers_1.9.0 latticeExtra_0.6-29    
## [67] png_0.1-7               vctrs_0.3.8             rmdformats_1.0.3       
## [70] gtable_0.3.0            purrr_0.3.4             assertthat_0.2.1       
## [73] xfun_0.28               lwgeom_0.2-8            e1071_1.7-9            
## [76] ragg_1.2.1              class_7.3-19            tibble_3.1.6           
## [79] ellipsis_0.3.2

Tutorial Part II — R goes spatial